*******		REPLICATION FILES
*******		Climate Variability and Irregular Migration to the European Union 
*******		Global Environmental Change 
*******		Fabien Cottier and Idean Salehyan
*******		Replication: Sensitivity analysis (S9) temperature and precipitation anomalies (Tables A.30-32, Figures A.22-24)
*******		This version: April 2021

* note: stata packages required for running models
* - coefplot
* - crossfold
* - plottig


clear all

set more off
set scheme plottig
cd "path/to/replication/directory"

* load data
use "data_quarter.dta", clear 

* load stata risk scripts
qui do scripts/script_rr_weather_082019

* duplicates data
duplicates list cowc year quarter
duplicates tag cowc year quarter, gen(_dupl)
duplicates drop cowc year quarter, force

* set time series
sort cowc year quarter
gen quarter_int = real(regexs(1)) if regexm(quarter,"([0-9]+)")
gen quarterS=(year-2005)*4+quarter_int
order cowc nationality year quarter quarter_int quarterS
tsset cowc quarterS
 
* gen additional variables
encode continent, gen(contN)
recode contN (5=.)


* generate lag quartely migration variables
by cowc: gen nmigrq_exbalk_1tl = nmigrq_exbalk[_n-1]
by cowc: gen nmigrq_exbalk_2tl = nmigrq_exbalk[_n-2]
by cowc: gen nmigrq_exbalk_3tl = nmigrq_exbalk[_n-3]
by cowc: gen nmigrq_exbalk_4tl = nmigrq_exbalk[_n-4]
by cowc: gen nmigrq_exbalk_ln_1tl = nmigrq_exbalk_ln[_n-1]
by cowc: gen nmigrq_exbalk_ln_2tl = nmigrq_exbalk_ln[_n-2]
by cowc: gen nmigrq_exbalk_ln_3tl = nmigrq_exbalk_ln[_n-3]
by cowc: gen nmigrq_exbalk_ln_4tl = nmigrq_exbalk_ln[_n-4]


* gen dummies variables for sample splitting: low agriculture / high agriculture
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* i.year if migr100_exbalk==1, ///
 fe vce(cluster cowc)

sum agriLab if e(sample) & year==2010, d
scalar agriMed=r(p50)
gen agri_H=0 if e(sample) & year==2010
recode agri_H (0=1) if agriLab-agriMed > 0 & year==2010
by cowc: replace agri_H = agri_H[21] if missing(agri_H)

* generate dummies variable for extreme weather (cut-off 10 percentile / 90 percentile)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* wmean_speiy /// 
 i.year i.quarter_int if migr100_exbalk==1 & year<=2015, ///
 fe vce(cluster cowc)

centile (stdw10ma_temp) if e(sample), centile (10 90)
gen low_tempy = 0 if stdw10ma_tempy !=.
recode low_tempy (0=1) if stdw10ma_tempy <= r(c_1)
centile (stdw10ma_temp) if e(sample), centile (10 90)
gen high_tempy = 0 if stdw10ma_tempy !=.
recode high_tempy (0=1) if stdw10ma_tempy >= r(c_2) 

centile (stdw10ma_precip) if e(sample), centile (10 90)
gen low_precipy = 0 if stdw10ma_precipy !=.
recode low_precipy (0=1) if stdw10ma_precipy <= r(c_1)
centile (stdw10ma_precip) if e(sample), centile (10 90)
gen high_precipy = 0 if stdw10ma_precipy !=.
recode high_precipy (0=1) if stdw10ma_precipy >= r(c_2)

*******		  Macros

global CLIM_Y0 stdw10ma_tempy stdw10ma_precipy 

global CLIM_SQ_Y0 c.stdw10ma_temp##c.stdw10ma_tempy c.stdw10ma_precip##c.stdw10ma_precipy 
 
global CLIM_Y0Y2 L(0 4 8).stdw10ma_tempy L(0 4 8).stdw10ma_precipy 

global CLIM_SQ_Y0Y2 c.stdw10ma_temp##c.stdw10ma_tempy cL4.stdw10ma_temp##cL4.stdw10ma_tempy cL8.stdw10ma_temp##cL8.stdw10ma_tempy /// 
  c.stdw10ma_precip##c.stdw10ma_precipy cL4.stdw10ma_precip##cL4.stdw10ma_precipy cL8.stdw10ma_precip##cL8.stdw10ma_preci

*******		 Table A.30

* NULL Model // No SPEI term
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0 /// 
	i.year i.quarter_int if migr100_exbalk==1 & year<=2015, ///
	fe vce(cluster cowc)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_*  /// 
	i.year i.quarter_int if e(sample), ///
	fe vce(cluster cowc)
est sto tabA30_m0
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse null model     : " CVrmse
* list countries included in sample 
est resto tabA30_m0
unique cowc if e(sample)
tabulate nationality contN if e(sample)


* Model 1
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0 /// 
	i.year i.quarter_int if migr100_exbalk==1 & year<=2015, ///
	fe vce(cluster cowc)
est sto tabA30_m1
* F test spei parameters & aic
test stdw10ma_tempy stdw10ma_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 2
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_SQ_Y0 ///
 i.year i.quarter_int if migr100_exbalk==1 & year<=2015, ///
 fe vce(cluster cowc)
est sto tabA30_m2
* F test spei parameters & aic
test stdw10ma_tempy c.stdw10ma_temp#c.stdw10ma_tempy stdw10ma_precipy c.stdw10ma_precip#c.stdw10ma_precipy 
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 3
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0Y2 /// 
	i.year i.quarter_int if migr100_exbalk==1 & year<=2015, ///
	fe vce(cluster cowc)
est sto tabA30_m3
* F test spei parameters & aic
test L0.stdw10ma_tempy L4.stdw10ma_tempy L8.stdw10ma_tempy ///
L0.stdw10ma_precipy L4.stdw10ma_precipy L8.stdw10ma_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 4
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_SQ_Y0Y2 ///
	i.year i.quarter_int if migr100_exbalk==1 & year<=2015, ///
	fe vce(cluster cowc)
est sto tabA30_m4
* F test spei parameters & aic
test L0.stdw10ma_tempy L4.stdw10ma_tempy L8.stdw10ma_tempy ///
c.stdw10ma_temp#c.stdw10ma_tempy cL4.stdw10ma_temp#cL4.stdw10ma_tempy cL8.stdw10ma_temp#cL8.stdw10ma_tempy ///
L0.stdw10ma_precipy L4.stdw10ma_precipy L8.stdw10ma_precipy ///
c.stdw10ma_precip#c.stdw10ma_precipy cL4.stdw10ma_precip#cL4.stdw10ma_precipy cL8.stdw10ma_precip#cL8.stdw10ma_precipy 
estat ic 
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse




*******		 Table A.31


* NULL Model // agrarian countries
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0 ///
	i.year i.quarter_int if migr100_exbalk==1 & agri_H==1 & year<=2015, ///
	fe vce(cluster cowc)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* ///
	i.year i.quarter_int if e(sample), ///
	fe vce(cluster cowc)
est sto tabA31_m0H
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse null model     : " CVrmse
* list of countries
est resto tabA31_m0H
unique cowc if e(sample)
tabulate nationality contN if e(sample)

* NULL Model // non-agrarian countries
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0 ///
	i.year i.quarter_int if migr100_exbalk==1 & agri_H==0 & year<=2015, ///
	fe vce(cluster cowc)
qui xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* ///
	i.year i.quarter_int if e(sample), ///
	fe vce(cluster cowc)
est sto tabA31_m0L
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse null model     : " CVrmse
global CVrmsenull_LA=round(CVrmse,0.001)
* list of countries
est resto tabA31_m0L
unique cowc if e(sample)
tabulate nationality contN if e(sample)


* Model 5 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0 ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==1 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m5
* F test spei parameters & aic
test stdw10ma_tempy stdw10ma_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 6 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==0 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m6
* F test spei parameters & aic
test stdw10ma_tempy stdw10ma_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 7 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_SQ_Y0 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==1 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m7
* F test spei parameters & aic
test stdw10ma_tempy c.stdw10ma_temp#c.stdw10ma_tempy stdw10ma_precipy c.stdw10ma_precip#c.stdw10ma_precipy 
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 8 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_SQ_Y0 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==0 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m8
* F test spei parameters & aic
test stdw10ma_tempy c.stdw10ma_temp#c.stdw10ma_tempy stdw10ma_precipy c.stdw10ma_precip#c.stdw10ma_precipy 
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 9 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0Y2 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==1 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m9
unique cowc if e(sample)
* F test spei parameters & aic
test L0.stdw10ma_tempy L4.stdw10ma_tempy L8.stdw10ma_tempy ///
L0.stdw10ma_precipy L4.stdw10ma_precipy L8.stdw10ma_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 10 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_Y0Y2 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==0 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m10
unique cowc if e(sample)
* F test spei parameters & aic
test L0.stdw10ma_tempy L4.stdw10ma_tempy L8.stdw10ma_tempy ///
L0.stdw10ma_precipy L4.stdw10ma_precipy L8.stdw10ma_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 11 // agrarian countriese
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_SQ_Y0Y2 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==1 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m11
* F test spei parameters & aic
test L0.stdw10ma_tempy L4.stdw10ma_tempy L8.stdw10ma_tempy ///
c.stdw10ma_temp#c.stdw10ma_tempy cL4.stdw10ma_temp#cL4.stdw10ma_tempy cL8.stdw10ma_temp#cL8.stdw10ma_tempy ///
L0.stdw10ma_precipy L4.stdw10ma_precipy L8.stdw10ma_precipy ///
c.stdw10ma_precip#c.stdw10ma_precipy cL4.stdw10ma_precip#cL4.stdw10ma_precipy cL8.stdw10ma_precip#cL8.stdw10ma_precip
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 12 // non-agrarian countries
* low reliance on agriculture
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* $CLIM_SQ_Y0Y2 ///
  i.year i.quarter_int  if migr100_exbalk==1 & agri_H==0 & year<=2015, ///
  fe vce(cluster cowc)
est sto tabA31_m12
* F test spei parameters & aic
test L0.stdw10ma_tempy L4.stdw10ma_tempy L8.stdw10ma_tempy ///
c.stdw10ma_temp#c.stdw10ma_tempy cL4.stdw10ma_temp#cL4.stdw10ma_tempy cL8.stdw10ma_temp#cL8.stdw10ma_tempy ///
L0.stdw10ma_precipy L4.stdw10ma_precipy L8.stdw10ma_precipy ///
c.stdw10ma_precip#c.stdw10ma_precipy cL4.stdw10ma_precip#cL4.stdw10ma_precipy cL8.stdw10ma_precip#cL8.stdw10ma_precip
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse



*******		 Table A.32



* Model 13 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* L(0).low_tempy L(0).high_tempy L(0).low_precipy L(0).high_precipy ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==1 & year<=2015, fe vce(cluster cowc)
est sto tabA32_m13
* F test spei parameters & aic
test low_tempy high_tempy low_precipy high_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 14 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* L(0).low_tempy L(0).high_tempy L(0).low_precipy L(0).high_precipy ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==0 & year<=2015, fe vce(cluster cowc)
est sto tabA32_m14
* F test spei parameters & aic
test low_tempy high_tempy low_precipy high_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse


* Model 15 // agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* L(0 4 8).low_tempy L(0 4 8).high_tempy L(0 4 8).low_precipy L(0 4 8).high_precipy ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==1 & year<=2015, fe vce(cluster cowc)
est sto tabA32_m15
* F test spei parameters & aic
test L0.low_tempy L4.low_tempy L8.low_tempy L0.high_tempy L4.high_tempy L8.high_tempy ///
L0.low_precipy L4.low_precipy L8.low_precipy L0.high_precipy L4.high_precipy L8.high_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse

* Model 16 // non-agrarian countries
xtreg nmigrq_exbalk_ln nmigrq_exbalk_ln_* L(0 4 8).low_tempy L(0 4 8).high_tempy L(0 4 8).low_precipy L(0 4 8).high_precipy ///
  i.year i.quarter_int if migr100_exbalk==1 & agri_H==0 & year<=2015, fe vce(cluster cowc)
est sto tabA32_m16
* F test spei parameters & aic
test L0.low_tempy L4.low_tempy L8.low_tempy L0.high_tempy L4.high_tempy L8.high_tempy ///
L0.low_precipy L4.low_precipy L8.low_precipy L0.high_precipy L4.high_precipy L8.high_precipy
estat ic
* cross validation / RMSE
local reg_call=e(cmdline)
qui crossfold `reg_call'
mata : st_numscalar("CVrmse", sqrt(mean((st_matrix("r(est)")):^2)))
di as text   "Avg rmse     : " CVrmse



*******		 Relative risk figures




* Figure A.22 (Model 2, Table A.30)

qui rr_est_clim tabA30_m2 "quadratic"


* Figure A.23 (Model 7 & 8, Table A.32) // Temperature anomalies

qui rr_est_clim_Sp tabA31_m7 tabA31_m8 "temp" "quadratic"


* Figure A.24 (Model 7 & 8, Table A.32) // Precipitation anomalies

qui rr_est_clim_Sp tabA31_m7 tabA31_m8 "precip" "quadratic"



*******		 Correlation SPEI and temperature / precipitation anomalies 


est resto tabA30_m1

* temperature
cor wmean_speiy stdw10ma_tempy if e(sample)

* precipitation
cor wmean_speiy stdw10ma_precipy if e(sample)

